NUSmodAn

Introduction

• Notes:
  1. AY2017/2018 Semester 2 and AY2018/2019 Semester 2 bidding data not available.
  2. The bidding statistics are highly non-normal due to being bounded by zero (they cannot make negative bids or have negative bidders). May consider using zero-inflated or poisson regression if considering these statistics as dependent variables.

Phase 1: Setting Up Environment, Packages And Loading Data

Packages And Options

library(semTools)
library(ggplot2)
library(rjson)
library(stringr)
library(DT)
library(psych)
library(corrplot)
library(tidyverse)
library(forcats)
library(lme4)
library(shiny)
library(semTools)
library(plotly)
options(width = 999)
knitr::opts_chunk$set(dpi = 300, out.width = "50%", eval = TRUE)

Load Bidding Data

• Extract data from nusmods at https://api.nusmods.com/.
• After loading the data in .JSON format, convert to a dataframe.

myjson <- fromJSON(file = url("https://api.nusmods.com/corsBiddingStatsRaw.json")) # read data directly from URL
myBid <- data.frame() # create empty dataframe which will act as a container to be populated with data
for(r in 1:length(myjson)) # for each element in the myjson list, append it to myBid
{
  if(myjson[[r]]$Semester == 1 | myjson[[r]]$Semester == 2) # if semester 1 or 2
  {
    myBid <- rbind(myBid, myjson[[r]]) # add to dataframe
  }
  myjson[[r]] <- NA # free up some RAM
}

saveRDS(myBid, file = "myBid.RDS") # save to directory

Load myBid.RDS

• Downloading the data from the API using the code above takes a substantial amount of time.
• I saved the downloaded data in myBid.RDS and load it directly from my local folder while I worked on the project.

myBid <- readRDS("myBid.RDS")

Load Module Information

• Module information was scattered across different folders.
• Used a loop to repeat the process of downloading and converting to dataframe across the different folders accessed by the different URLs.
  • The same concept was used to consolidate information about the Module Titles.

myModInfo <- data.frame() # create empty dataframe which will act as a container to be populated with data
for(year in c(2011:2018)) # looping through each year
{
  for(semester in c(1,2))
  {
    # create the url where data is to be extracted from
    myurl <- paste0("https://api.nusmods.com/", year, "-", year + 1, "/", semester, "/moduleTimetableDeltaRaw.json")
    myjson <- fromJSON(file = url(myurl))
    for(r in 1:length(myjson)) # for each element in the myjson list, append it to myModInfo
    {
      if(isTRUE(str_detect(myjson[[r]]$ModuleCode, "^PL"))) # only keep info if module code begins with PL
      {
        if(myjson[[r]]$Semester == 1 | myjson[[r]]$Semester == 2) # only get semester 1 and 2 information
        {
          myModInfo <- rbind(myModInfo, myjson[[r]]) # add to dataframe
        }
      }
      myjson[[r]] <- NA # replace the element with NA to free up some rAM
    }
    cat(year, "Semester", semester, "Done!") # progress tracker
  }
}

myTitles <- data.frame() # create empty dataframe which will act as a container to be populated with data
for(year in c(2014:2018)) # looping through each year
{
    myurl <- paste0("https://api.nusmods.com/", year, "-", year + 1, "/moduleList.json") # create the url where data is to be extracted from
    myjson <- fromJSON(file = url(myurl))
    for(r in 1:length(myjson)) # for each element in the myjson list, append it to myModInfo
    {
      if(isTRUE(str_detect(myjson[[r]]$ModuleCode, "^PL"))) # only keep info if module code begins with PL
      {
        if(paste0(myjson[[r]]$Semester, collapse = "|") == "1"|
           paste0(myjson[[r]]$Semester, collapse = "|") == "2"|
           paste0(myjson[[r]]$Semester, collapse = "|") == "1|2") # only keep information from semester 1 and 2
        {
          myTitles <- rbind(myTitles, as.data.frame(myjson[[r]])) # add to dataframe
        }
      }
      myjson[[r]] <- NA # free RAM
    }
}

myModInfo <- myTitles %>% # add titles information to myModInfo
  select(ModuleCode, ModuleTitle) %>% # select these two columns
  filter(ModuleTitle != "Lab in Applied Psychology") %>%
  distinct() %>% # remove duplicates
  right_join(myModInfo, by = "ModuleCode") # left = myTitles, right = myModInfo

saveRDS(myModInfo, file = "myModInfo.RDS") # save to directory

Load myModInfo.RDS

• Downloading the data from the API using the code above takes a substantial amount of time.
• I saved the downloaded data in myModInfo.RDS and load the data directly while I worked on the project.

myModInfo <- readRDS("myModInfo.RDS")

Phase 2: Filter, Transform And Merge

Module Information

• Filter information from the dataframe myModInfo.
  • Removing non-Psychology modules.
  • Removing modules without module titles, these are modules that appeared before AY2014/2015 and never resurfaced afterwards.
  • Removing information about tutorials.

Filter

myModInfo <- myModInfo %>%
  select(-LastModified, -LastModified_js, -isDelete) %>% # remove these columns
  filter(str_detect(ModuleCode, "^PL")) %>% # removing non-Psychology modules
  filter(!is.na(ModuleTitle)) %>% # removing modules without module titles #PL3285, PL4220, PL4217
  filter(LessonType != "TUTORIAL") %>% # removing information about tutorials
  select(AcadYear, Semester, ModuleCode, ModuleTitle, DayText, StartTime, Semester) %>% # select these columns
  distinct() # remove duplicates

head(myModInfo) # peek

Bidding Information myBid

• Filter information from the dataframe myBid.
  • Removing non-Psychology modules, including Roots and Wings (prefixed with PLS-) and Psychology for non-Psychology students (prefixed with PLB-).
  • Removing information from quotas that are reserved and not available for bidding.
  • Removing information from modules with more than one lecture/seminar session.
  • Removing bidding information from non-psychology students.

Filter

myBid <- myBid %>%
  filter(str_detect(ModuleCode, "^PL")) %>% # removing non-Psychology modules
  filter(!str_detect(ModuleCode, "PLS|PLB")) %>% # remove PLS and PLB modules
  filter(!str_detect(StudentAcctType, "Reserved")) %>% # remove reserved rounds
  filter(!str_detect(StudentAcctType, "[G]")) %>% # remove bidding information from non-psychology students
  filter(!str_detect(paste0(unique(myBid$ModuleCode[grep("2",myBid$Group)]), collapse = "|"), ModuleCode)) %>% # remove modules that have more than one lecture
  select(-Faculty, -Group) # remove these columns

head(myBid) # peek

Merge

• Combine the information of myModInfo and myBid.

# modules that do not appear in both dataframes are dropped
mydata <- inner_join(myBid, 
                     myModInfo,
                     by = c("ModuleCode", "AcadYear", "Semester"))

Phase 3: Data Wrangling

• The variables available in the original data are useful but they are too specific to interpret meaningfully.
• This section creates new variables based on the original data and allow us to better discern any trend in the data.
• Also includes additional wrangling and manipulations to ease the plotting of graphs and analysis later.

Coercing Columns

# transform these columns to numeric
for(r in c("Quota", "Bidders", "LowestBid", "LowestSuccessfulBid", "HighestBid", "StartTime"))
{
  mydata[,grep(r, names(mydata))] <- as.numeric(mydata[,grep(r, names(mydata))])
}
# transform these columns to factors
for(r in c("AcadYear", "Semester", "ModuleCode", "Round", "StudentAcctType", "DayText", "StudentAcctType", "ModuleTitle"))
{
  mydata[,grep(r, names(mydata))] <- factor(mydata[,grep(r, names(mydata))])
}

Rearranging DayText Levels

mydata$DayText <- factor(mydata$DayText,
                         levels = c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday"))

Shortening StudentAcctType Levels

mydata <- mutate(mydata, StudentAcctType = fct_recode(StudentAcctType, 
                                                      "New[P]" = "New Students [P]",
                                                      "NUS[P]" = "NUS Students [P]",
                                                      "Return[P]" = "Returning Students [P]",
                                                      "ReturnNew[P]" = "Returning Students and New Students [P]"))

New Variable: Level

# create new variable that indicates the level of the module, based on their module code
mydata$Level <- factor(ifelse(str_detect(mydata$ModuleCode, "1[0-9][0-9][0-9]"), "Level 1",
                      ifelse(str_detect(mydata$ModuleCode, "2[0-9][0-9][0-9]"), "Level 2",
                             ifelse(str_detect(mydata$ModuleCode, "3[0-9][0-9][0-9]"), "Level 3",
                                    ifelse(str_detect(mydata$ModuleCode, "4[0-9][0-9][0-9]"), "Level 4", 
                                           "Graduate Module")))))

New Variable: BidPerQuota

# BidPerQuota = Bidders/Quota
mydata$BidPerQuota <- with(mydata, Bidders/Quota)

New Variable: DayPeriod

mydata$Period <- factor(ifelse(mydata$StartTime < 1200, "Morning", ">=Afternoon"),
                           levels = c("Morning", ">=Afternoon"))

Vector Of Column Names

# create vector of the column names which are factors
facnames <- mydata %>% select_if(is.factor) %>% names()
# facnames without ModuleCode and StudentAcctType
facnames.mod <- facnames[-grep("ModuleCode|ModuleTitle", facnames)]
# create vector ofthe column names which are numeric
numnames <- mydata %>% select_if(is.numeric) %>% names()
# numnames without StartTime
numnames.time <- names(select_if(mydata, is.numeric))[-grep("StartTime", numnames)]

Phase 4: Data Diagnostics

• Plot univariate histograms and bivariate plots using loops for almost every combination of variables.
• The graphs from this section are predominantly for diagnostics rather than exploration, what I mean is that the graphs from this section would make little sense if one tried to draw insights from them. This is because they are aggregated across all other variables.
  • For example: The mean of Bidders is calculated across all academic years, all bidding rounds, all modules…
• What I am looking out for in this section are odd patterns, like zeroes in places where they shouldn’t be, missing data, highly non-normal data, variables with outliers, etc…

Univariate Descriptive Statistics

str(mydata)

## 'data.frame':    1621 obs. of  16 variables:
##  $ AcadYear           : Factor w/ 8 levels "2011/2012","2012/2013",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Semester           : Factor w/ 2 levels "1","2": 2 2 2 2 2 2 2 2 2 2 ...
##  $ Round              : Factor w/ 7 levels "1A","1B","1C",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ ModuleCode         : Factor w/ 76 levels "PL2131","PL2132",..: 1 1 2 2 3 3 4 4 5 5 ...
##  $ Quota              : num  5 12 35 35 28 50 25 22 25 30 ...
##  $ Bidders            : num  3 42 8 3 7 2 8 5 3 3 ...
##  $ LowestBid          : num  1 205 1 1 1 1 1 1 1 1 ...
##  $ LowestSuccessfulBid: num  1 977 1 1 1 1 1 1 1 1 ...
##  $ HighestBid         : num  368 1255 500 250 1200 ...
##  $ StudentAcctType    : Factor w/ 4 levels "New[P]","NUS[P]",..: 3 1 3 1 3 1 3 1 3 1 ...
##  $ ModuleTitle        : Factor w/ 74 levels "Abnormal Psychology",..: 64 64 65 65 7 7 12 12 16 16 ...
##  $ DayText            : Factor w/ 5 levels "Monday","Tuesday",..: 3 3 2 2 2 2 3 3 1 1 ...
##  $ StartTime          : num  1600 1600 800 800 1200 1200 1400 1400 1400 1400 ...
##  $ Level              : Factor w/ 3 levels "Level 2","Level 3",..: 1 1 1 1 2 2 2 2 2 2 ...
##  $ BidPerQuota        : num  0.6 3.5 0.2286 0.0857 0.25 ...
##  $ Period             : Factor w/ 2 levels "Morning",">=Afternoon": 2 2 1 1 2 2 2 2 2 2 ...

select(describe(mydata), n, mean, sd, median, min, max, skew, kurtosis)

Univariate Histograms/Boxplots

• Plots to illustrate the frequency/distribution of variables.

Categorical Variables

# plot the categorical variables
for(r in facnames.mod)
{
  plot(
    ggplot(data = mydata, aes_string(x = r, fill = r)) + 
      geom_histogram(stat = "count") + 
      ylab("Count") +
      ggtitle(paste0("Count of ", r)) +
      theme_classic() + 
      theme(axis.text.x = element_text(angle = 90),
            axis.title.x = element_blank(),
            legend.position = "none")
  )
}

Continuous Variables

# plot the continuous variables
for(r in numnames)
{
  plot(
    ggplot(data = mydata, aes_string(x = r, fill = r)) + 
      geom_histogram(fill = "violetred", alpha = 0.5, bins = 50) + 
      ylab("Histogram") +
      ggtitle(paste0(r)) +
      theme_classic() + 
      theme(axis.text.x = element_text(angle = 90),
            axis.title.x = element_text())
  )
}

Bivariate Plots

• Plots to illustrate pairwise relationships amongst variables.

Categorical-Categorical

for(r in 1:length(facnames.mod)) # loop across all factors
{
  for(i in 1:length(facnames.mod)) # inner loop
  {
    if(i == r | i < r) 
    { # dont do anything if they are the same or the graph has been made before
    } else {
      tempform <- paste0("~ ", facnames.mod[r], " + ", facnames.mod[i])  # create formula for xtabs
      # temp is a dataframe that is only going to exist in this section and overwritten with each loop
      temp <- as.data.frame(xtabs(eval(parse(text = tempform)),
                                  data = mydata,
                                  subset = NULL))
      plot(
        ggplot(data = temp, aes_string(x = facnames.mod[r], y = facnames.mod[i], fill = "Freq", label = "Freq")) +
          geom_tile() + 
          geom_text() + 
          scale_fill_gradient(low = "white", high = "violetred") + 
          theme_minimal() + 
          theme(axis.text.x = element_text(angle = 90),
                legend.position = "none")
      )
    }
  }
}

Continuous-Continuous

for(r in 1:length(numnames)) # loop across all numeric columns
{
  for(i in 1:length(numnames)) # inner loop
  {
    if(i == r | i < r)
    { # dont do anything if they are the same or the graph has been made before
    } else {
      # create formulas for lm()
      tempform.std <- paste0("scale(", numnames[i],")", " ~ ", "scale(", numnames[r], ")") # standardized
      tempform <- paste0(numnames[i], " ~ ", numnames[r]) # unstandardized
      # regress to get best fit line
      stdreg <- lm(eval(parse(text = tempform.std)),
                   data = mydata) # standardized
      reg <- lm(eval(parse(text = tempform)),
                data = mydata) # unstandardized
      
      plot(
        ggplot(data = mydata, aes_string(x = numnames[r], y = numnames[i])) +
          geom_point(color = "violetred", size = 2, alpha = 0.3) +
          theme_classic() + 
          geom_abline(slope = reg$coefficients[2], intercept = reg$coefficients[1], lty = "dashed") + 
          geom_label(aes(x = Inf, y = Inf, label = paste0("Standardized Regression Coefficient = ",
                                                          round(stdreg$coefficients[2],3)),
                         hjust = 1, vjust = 2)) + 
          theme(axis.text.x = element_text(angle = 90))
      )
    }
  }
}

Continuous-Categorical

for(r in facnames.mod) # loop across all factor columns
{
  for(i in numnames) # inner loop across all numeric columns
  {
    plot(
      ggplot(data = mydata, aes_string(x = r, y = i, fill = r)) + 
        geom_boxplot() + 
        theme_classic() + 
        theme(legend.position = "none",
              axis.text.x = element_text(angle = 90))
    )
  }
}

Correlation Matrix

corrplot.mixed(cor(mydata[,grep(paste0(numnames.time, collapse = "|"), names(mydata))]), # only get the numeric columns
               upper = "color",
               tl.pos = "lt",
               tl.cex = 0.5,
               cl.cex = 0.5)

Phase 5: Answering Questions

Is it true that it is easier to bid for a module in the morning?

Do bids become higher as the rounds get later?

Most/Least Popular Level 4 Module

What are the most and least popular modules?

Lowest Successful Bid

ggplotly(
  mydata %>%
  filter(Level ==  "Level 4") %>%
  filter(Round == "1A") %>%
  group_by(ModuleCode, ModuleTitle) %>%
  summarize(LSB.avg.sem.years = mean(LowestSuccessfulBid)) %>%
  ungroup() %>%
  mutate(ModuleCode = fct_reorder(ModuleCode, LSB.avg.sem.years)) %>%
  ggplot(mapping = aes(x = ModuleCode, y = LSB.avg.sem.years, label = ModuleTitle, fill = LSB.avg.sem.years)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  theme_classic() +
  theme(legend.position = "none")
, tooltip = c("x", "label", "y")
, height = 600, width = 400
)

BPQ

ggplotly(
  mydata %>%
  filter(Level ==  "Level 4") %>%
  filter(Round == "1A") %>%
  group_by(ModuleCode, ModuleTitle) %>%
  summarize(BPQ.avg.sem.years = mean(BidPerQuota)) %>%
  ungroup() %>%
  mutate(ModuleCode = fct_reorder(ModuleCode, BPQ.avg.sem.years)) %>%
  ggplot(mapping = aes(x = ModuleCode, y = BPQ.avg.sem.years, label = ModuleTitle, fill = BPQ.avg.sem.years)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  theme_classic() +
  theme(legend.position = "none")
, tooltip = c("x", "label", "y")
, height = 600, width = 400
)

Quota

ggplotly(
  mydata %>%
  filter(Level ==  "Level 4") %>%
  filter(Round == "1A") %>%
  group_by(ModuleCode, ModuleTitle) %>%
  summarize(quota.avg.sem.years = mean(Quota)) %>%
  ungroup() %>%
  mutate(ModuleCode = fct_reorder(ModuleCode, quota.avg.sem.years)) %>%
  ggplot(mapping = aes(x = ModuleCode, y = quota.avg.sem.years, label = ModuleTitle, fill = quota.avg.sem.years)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  theme_classic() +
  theme(legend.position = "none")
, tooltip = c("x", "label", "y")
, height = 600, width = 400
)

Do less people bid for a module if the lecture begins in the morning (before 12pm)?

Lets look at each module and compare the average number of bidders, bidders per quota and lowest successful bids when the lecture begins in and after the morning.

By Module

for(r in c("meanBidders", "meanBpQ", "meanLSB"))
{
  plot(mydata %>%
    group_by(ModuleCode, ModuleTitle, Period) %>%
    summarise(meanBidders = mean(Bidders), meanBpQ = mean(BidPerQuota), meanLSB = mean(LowestSuccessfulBid),
              sdBidders = sd(Bidders), sdBpQ = sd(BidPerQuota), sdLSB = mean(LowestSuccessfulBid)) %>%
    ggplot(aes_string(x = "Period", y = r, fill = "Period")) +
    geom_bar(stat = "identity") +
    theme_classic() +  
    theme(axis.text.x = element_blank(), 
          legend.position = "top",
          strip.background = element_rect(fill = "grey30", linetype = "blank"),
          strip.text = element_text(color = "white", size = 9)) + 
    facet_wrap(~ ModuleCode:ModuleTitle, labeller = label_wrap_gen(width = 25)) +
    ggtitle(r))
}

By Level

for(r in c("meanBidders", "meanBpQ", "meanLSB"))
{
  plot(mydata %>%
    group_by(Level, Period) %>%
    summarise(meanBidders = mean(Bidders), meanBpQ = mean(BidPerQuota), meanLSB = mean(LowestSuccessfulBid),
              sdBidders = sd(Bidders), sdBpQ = sd(BidPerQuota), sdLSB = mean(LowestSuccessfulBid)) %>%
    ggplot(aes_string(x = "Period", y = r, fill = "Period")) +
    geom_bar(stat = "identity") +
    theme_classic() +  
    theme(axis.text.x = element_blank(), 
          legend.position = "top",
          strip.background = element_rect(fill = "grey30", color = "black"),
          strip.text = element_text(color = "white", size = 12)) + 
    facet_wrap(~ Level) +
    ggtitle(r))
}

Bonus: Multilevel Modeling

Removal Of Modules With…

Less Than 3 Instances

mymlm <- mydata %>%
  group_by(ModuleCode) %>%
  filter(n() >= 3)

With No Variability In StartTime

# # remove those with no variability in StartTime
mymlm <- mydata %>%
  group_by(ModuleCode) %>%
  mutate(sdBpQ = sd(BidPerQuota), sdStartTime = sd(StartTime)) %>%
  filter(sdStartTime > 0) %>%
  arrange(sdBpQ)

Hours From 12am Hrs12 And Group-Mean Centered Hrs12.

mymlm <- mymlm %>% 
  mutate(Hrs12 = StartTime/100) %>%
  group_by(ModuleCode) %>%
  arrange()

Peek Data

head(mymlm)

Models

Baseline Model

fit0 <- lmer(BidPerQuota ~ 1 + (1 | ModuleCode),
     data = mymlm)
VarCorr(fit0)

(0.3672^2)/((0.3672^2) + (0.67259^2))

With Hrs12.gmc

fit1 <- lmer(avgLSB ~ Level*Hrs12.gmc + (1 | ModuleCode),
     data = mymlm)
summary(fit1)
car::Anova(fit1)

Do results from previous rounds…

head(mydata)

a <- mydata %>%
  pivot_wider(names_from = Round,
         values_from = BidPerQuota
         )

are there dumping